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Abstract 



Using a new, second-order accurate numerical method we present dynamical simula- 
tions of oblique MHD cosmic ray (CR) modified plane shock evolution using the two-fluid 
model for diffusive particle acceleration. The numerical shocks evolve to published analyti- 
cal steady state properties. In order to probe the dynamical role of magnetic fields we have 
explored for these time asymptotic states the parameter space of upstream fast mode Mach 
number, Mf, and plasma f3, compiling the results into maps of dynamical steady state CR 
acceleration efficiency, e c . These maps, along with additional numerical experiments, show 
that e c is reduced through the action of compressive work on tangential magnetic fields 
in CR-MHD shocks. Thus e c in low (3, moderate Mf shocks tends to be smaller in quasi 
perpendicular shocks than it would be high (3 shocks of the same My. This result supports 
earlier conclusions that strong, oblique magnetic fields inhibit diffusive shock acceleration. 
For quasi parallel shocks with (3 < 1, on the other hand, e c seems to be increased at a given 
Mf when compared to high (3 shocks. The apparent contradiction to the first conclusion 
results, however, from the fact that for small (3 quasi parallel shocks, the fast mode Mach 
number is not a good measure of compression through the shock. That is better reflected 
in the sonic Mach number, which is greater. Acceleration efficiencies for high and low f3 
having comparable sonic Mach numbers are more similar. 



Time evolution of CR-MHD shocks is qualitatively similar to CR-gasdynamical shocks. 
However, several potentially interesting differences are apparent. We have run simulations 
using constant, and non-isotropic, obliquity dependent forms of the diffusion coefficient k. 
Comparison of the results show that while the final steady states achieved are the same 
in each case, the history of CR-MHD shocks can be strongly modified by variations in 
K and, therefore, in the acceleration timescale. Also, the coupling of CR and MHD in 
low f3, oblique shocks substantially influences the transient density spike that forms in 
strongly CR-modified shocks. We find that inside the density spike a MHD slow mode 
wave is generated that eventually steepens into a shock. A strong shear layer develops 
within the density spike, driven by MHD stresses. We conjecture that currents in the 
shear layer could, in non-planar flows, result in enhanced particle acceleration through 
drift acceleration. 
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1. Introduction 

In the past several years nonlinear theories of diffusive shock acceleration have demon- 
strated the potential importance of cosmic-ray (CR) dynamical feedback on the evolution 
of shock structures (see reviews by Blandford & Eichler 1987; Berezhko & Krymskii 1988; 
Jones & Ellison 1991). A range of methods has been successfully applied to show that CR 
pressures can become sufficient to modify shock structures substantially (e.g., Drury & 
Volk 1981; Achterberg et al. 1984; Ellison & Eichler 1984; Webb, Drury, and Volk 1986; 
Falle & Giddings 1987; Kang & Jones 1990, 1991). In some cases the modifications can 
even produce smooth CR dominated shocks with no entropy generating gas sub-shock. Nu- 
merical simulations using a variety of techniques have also demonstrated the importance of 
time-dependent effects in the determination of CR modified shock properties (e.g., Drury 
& Falle 1987; Falle & Giddings 1989; Kang & Jones 1991; Duffy 1992; Kang, Jones, & 
Ryu, 1992; Ryu, Kang & Jones 1993). 

Through all of this the presence of a large-scale magnetic field was generally understood 
to be an essential physical ingredient in CR shocks, since magnetic field fluctuations (small- 
scale Alfven waves, or "Alfvenic turbulence" to avoid confusion with large scale waves) 
mediate the CR propagation. But the magnetic fields are difficult to treat explicitly and 
fully, especially when those fields play a direct dynamical role. So, most studies of CR 
modified shocks have focused on pure gasdynamical models of the flows. That may be 
adequate if the fields are sufficiently weak, and/or the large-scale field is aligned parallel 
to the shock normal, so that sonic-mode shock dynamics applies. Except for modeling the 
spatial dependence of the CR diffusion coefficient, explicit treatment of the magnetic field 
may not always be generally necessary when the shock speed in the preshock (upstream) 
fluid, u s , and the preshock sound speed, a are large compared to the preshock Alfven speed, 
b. The second of these conditions is commonly expressed as p = a jb ^> 1. But, this 
condition is likely to be violated, perhaps strongly, in many environments where particle 
acceleration occurs (e.g., Slavin & Cox 1992). Even if the magnetic field and shock normal 
are parallel, treatment of CR shocks in flows with f3<l may depend in significant ways 
on /3, especially if the Alfvenic Mach number is not very large. One must be concerned, 
for example, with the finite momentum and energy carried by the CR-mediating Alfvenic 
turbulence (e.g., Volk, Drury & McKenzie 1984; Jones 1993a) 

Postponing for now discussion of the more complicated topics in general magnetohy- 
drodynamical (MHD) flows, our purpose in the present paper is to examine some straight- 
forward dynamical influences of oblique magnetic fields on plane CR shocks. There are 
several obvious ways in which oblique fields may be important. If the pressure in the tan- 
gential field is comparable to the gas pressure of the incoming flow, that changes the jump 
conditions to be satisfied through the shock and must, therefore, influence the transfer 
of energy and momentum to CR within the shock (e.g., Webb 1983). The anisotropic 
nature of MHD stresses in oblique shocks leads to the existence of fast, intermediate and 
slow wave mode families, each with distinct behaviors. That can complicate the shock 
properties and lead to the generation of distinctly MHD time dependent structures in a 
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shock evolving in response to CR acceleration. This increase in the degrees of freedom at 
the shock also augments the modes through which the shock structure may become dy- 
namically unstable (e.g., Chalov 1988; Zank et al. 1990). The CR diffusion coefficient may 
be anisotropic relative to the field, affecting the rates at which CR are accelerated. If the 
diffusion is described by standard kinetic theory models, the acceleration rate at oblique 
shocks could be much faster than at otherwise comparable parallel shocks (Jokipii 1987; 
Ostrowski 1988), although the resulting time-asymptotic test-particle form for the particle 
spectrum is apparently unchanged (Bell 1978), at least for nonrelativistic shocks (Kirk 
& Heavens 1989). This change in the acceleration rate should influence the dynamical 
properties of the shock at intermediate times, before a dynamical equilibrium is achieved, 
even if the field is weak and not directly involved dynamically. On the other hand, Baring, 
Ellison & Jones (1993) have shown that the efficiency of thermal particle injection into the 
CR population can be decreased by perpendicular magnetic field components. From these 
examples it is clear that as a step to understand the fundamental nature of diffusive shock 
acceleration in more realistic astrophysical settings, full time dependent MHD calculations 
are needed. 

Two-fluid models of CR transport along the lines developed by Drury & Volk (1981) 
and others are an efficient means to begin such explorations. These models treat transport 
of the CR through an energy conservation equation derived from the kinetic, diffusion- 
advection equation (e.g., Skilling 1975). The CR momentum distribution function itself 
is not followed in detail, except possibly through a priori models (e.g., Jones & Kang 
1993; Duffy et al. 1993). Two-fluid methods are meant to be used primarily for dynamical 
studies. They enable one to calculate economically the dynamical features of flows within 
the constraints imposed by the need to estimate a priori some closure parameters for the 
CR (see Drury 1983; Kang & Jones 1991; Kang 1993). In particular one must model 
the adiabatic index for the CR and the energy- weighted mean diffusion coefficient. To 
the extent that the solutions depend on those parameters (see, e.g., Achterberg et al. 
1984; Kang & Jones 1990) in ways that are not modelled, two-fluid models should not 
be used for quantitative prediction. On the other hand, the dependence of two-fluid 
solutions on these parameters for gasdynamic CR problems has been reasonably well 
studied (e.g., Achterberg et al. 1984; Kang & Jones 1990) and the basic validity of two-fluid 
dynamical solutions compared to more complete, diffusion-advection solutions for the CR 
momentum distribution demonstrated (e.g., Falle & Giddings 1987; Kang & Jones 1991; 
Ryu, Kang & Jones 1993). Some recent calculations even suggest that relatively simple 
models connecting the properties of the particle distribution to the closure parameters 
might provide adequate and internally self-consistent procedures for eliminating most 
limitations (Duffy et al. 1993). In any case, our purposes here are to explore the dynamical 
consequences of the magnetic field in these problems. Two-fluid methods would seem to 
be both adequate and appropriate for that aim. 

Webb (1983) first applied two-fluid methods to equilibrium, steady state CR-MHD 
modified shock structures for the case where the upstream gas is cold (its pressure can be 
ignored) and the magnetic field is oblique. Webb's analytic calculations demonstrated that, 
as for the CR modified gasdynamical flows, both sub-shock containing and smoothed, CR 
dominated solutions to the CR-MHD shock equations were possible. However, his solutions 
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suggested that shock CR acceleration was less effective when the upstream tangential 
component of the field was strong. Subsequently Kennel, Edmiston & Hada (1985) and 
Webb, Drury, and Volk 1986 (hereafter, WDV) calculated steady state two-fluid CR-MHD 
modified shock structures with finite upstream gas pressures. From their analysis, which 
included (3 = 5 and (3 = 1 cases and a range of fast mode Mach numbers, My, they 
concluded that the efficiency of CR acceleration was suppressed when the obliquity of the 
upstream magnetic field was high. They also found that the reduction in CR acceleration 
efficiency was most pronounced in the lower (3 ((3 = 1) flows of low fast mode Mach number 
(My < 6). For high (3 {(3 = 5) flows the CR acceleration efficiency was only weakly affected 
by changes in the upstream magnetic field angle regardless of the fast mode Mach number. 

As an extension of the previous work, we report in this paper the results of time 
dependent CR-MHD two-fluid simulations using a new numerical code. In an earlier 
work (Frank, Jones, & Ryu 1993) we presented tests of this code against Webb's cold 
gas analytical models as well as more general internal checks on our code's ability to 
evolve oblique shocks to self-consistent steady state CR-MHD structures. The simulations 
in that paper also confirmed Webb's cold gas solutions, presented preliminary results from 
more general CR-MHD flows and determined the convergence properties of our numerical 
method. In the present paper we explore the time evolution of finite gas pressure (Pg ^ 0) 
models of CR-MHD modified oblique shocks and their relation to the WDV steady state 
calculations. We will also explore in more detail questions raised in WDV concerning the 
role of both the obliquity of the field and the tangential field strength in determining the CR 
acceleration efficiency. More general time dependent effects unique to CR-MHD modified 
shocks will also be examined along with their potential importance for CR studies. As in 
WDV we will concentrate on those flows with high enough Alfvenic Mach number that 
switch-on and switch-off shocks do not occur. In the regime where those special shocks 
form, CR scattering centers (the Alfvenic turbulence) move at speeds comparable to the 
gas, requiring that the time dependent wave field dynamics be included in the calculations. 
For the present we restrict ourselves to flows where the scattering centers can be considered 
moving with the fluid. Jun, Clarke & Norman (1993) have also recently simulated CR- 
MHD shocks using two-fluid methods. Their results are complementary to ours, since they 
restricted their analysis to the special case of perpendicular (magnetosonic) shocks, so that 
they could not explore a number of questions address here. 

The remainder of the paper is organized as follows. Section 2 introduces the funda- 
mental equations governing CR modified MHD flows and discusses their application to 
our methods. An appendix at the end discusses these numerical methods in more detail. 
Section 3 compares the results of our models with those of WDV and explores issues con- 
cerning the dependence of CR shock acceleration on upstream magnetic field properties. 
In §4 we examine the general time dependent features of CR-MHD modified shocks as seen 
in our simulations, while we provide in §5 a discussion of these results in light of previous 
studies and their potential for application in higher dimensional flows. We present our 
conclusions in 56. 



2. Methods 
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We solve the equations of non-relativistic ideal MHD for one dimensional flow in 
Cartesian coordinates (e.g., Jeffrey 1968). As with gasdynamical models the conservation 
equations are modified to include momentum and energy source terms from CR feedback 
(e.g., Drury & Volk, 1981; Jones & Kang, 1990). The MHD equations are written in 
conservative, vector form as 

$ + ^ = S, (2.1) 



with 



and 



F = 



dt dx 
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pu z 
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(2.3) 



while the CR source term vector is 



S= - 



( o \ 

dPc/dx 





\u x dP c /dx + S e J 



(2.4) 



In two-fluid models the CR are themselves treated as a massless, diffusive fluid through a 
conservation equation for the CR energy, E c , derived by taking the second moment of the 
diffusion-advection equation (e.g., Achterberg 1982); namely, 



dE c d_ ( dEc 
dt dx 1 x c dx 



- P, 



du x 
dx 



+ S e 



(2.5) 



A more general treatment would replace equation [2.5] with the diffusion-advection equa- 
tion itself. Its moments would be calculated independently. In these relations 



P* = Pg + \(B 2 X + B 2 y + Bl) =P g + l -B\ 



(2.6) 



E = \p(u\ + u 2 y + u\) + -j— -— jj-P + \(B 2 X + B 2 y + Bl), 
P c = ( 7c - 1)E C , 



(2.7) 
(2.8) 
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and the magnetic field components are expressed in rationalized units 



B ^ 7S- (2 ' 9) 



so that the Alfven speed, b = B / \fp. In the expressions presented above the following 
definitions hold: p is the mass density; u x , B x and u y , u z , By, B z are the components of 
velocity and magnetic field parallel and perpendicular to shock front; P g , 7 5 and _P C , 7 C 
are the gas and CR pressures and adiabatic indices, and (k) is the energy-weighted, CR 
diffusion coefficient parallel to the shock normal (Drury 1983). (For simplicity, we will 
henceforth omit the brackets around k.) It is not necessary in two-fluid models to assume 
that either of 7 C or k, is a constant or that k derives from an energy independent diffusion 
process. As described earlier, they do depend on a priori models, however. The quantity 
S e is a term that allows direct energy transfer from the gas to the CR, such as through 
the injection of thermal particles into the CR population (e.g., Jones & Kang 1990). This 
is introduced for completeness, but for the present, we set S e = 0. We note in passing, 
however, that in the course of this investigation we carried out analogous calculations with 
S e 7^ 0. Behaviors of these shocks were not qualitatively different from those we will 
describe (cf. Kang & Jones 1990, 1991). In this discussion 7 5 = 5 /3, while the CR closure 
parameter 7 C , which depends upon the mix of non-relativistic (7 C — > 5 /3) and relativistic 
(7c — > 4 /3) particles in the CR population, will be treated as an input parameter. In 
general both 7 C and k are properties of the solution, so the need to specify their values a 
priori is the principal restriction of the two-fluid model. 



In the interest of simplicity we will not consider here flows with rotations of the shock- 
normal/magnetic-field plane between the initial upstream and downstream states, although 
the numerical code is quite capable of handling such features. Thus, without further 
loss of generality we can place the magnetic field in the X-Z plane, B = (B x , 0, B z ) = 
B(cos9 : 0, sin9). We will also select a reference frame so that there is no upstream 
tangential velocity; u y = u z = 0. All of the simulations discussed here are piston driven 
shock tubes. We establish the flows by projecting magnetized fluid with embedded CR 
in from the right boundary, using continuous boundary conditions and reflecting it off a 
wall (piston) at the left boundary. The tangential magnetic field at the left boundary 
is "mirrored" in a manner analogous to the gas density. Thus, there is no current sheet 
on the piston surface. In other words, the simulations can be considered as those of two 
identical, colliding clouds with magnetic field "continuous" across them. The CR pressure 
is given continuous boundary conditions on both ends of the grid. 

In general, spatial diffusion will differ in directions parallel and perpendicular to the 
magnetic field. Diffusive acceleration depends on the diffusion along the shock normal, 
since that determines the rate at which particles re-cross the shock. Thus, the appropriate 
diffusion coefficient in equation [2.5] will take the form (e.g., Drury 1983; Webb 1983), 

k = k\\cos 2 + K ± sm 2 6, (2.10) 

where the directions || and _L refer to the magnetic field direction. Note, of course, that 9 
is generally a changing function of space and time. In standard kinetic theory scattering 



models of particle diffusion the ratio ~ 1 + (<ut) , where w, and r are the particle 

gyro frequency and collision time, respectively (e.g., Jokipii 1987; Zank et al. 1990). We 
will adopt the form for k in equation [2.10] in our simulations. 

Previous studies of CR modified shocks have shown that the time required to evolve 
towards dynamical equilibrium scales with the so-called diffusion time, t^, (e.g., Drury & 
Falle 1986; Jones & Kang 1990), which in the present case is conveniently expressed as 

*d=4 = -. ( 2 - n ) 

u% u s 

where u s is the shock speed in the upstream fluid (see equation [3.7]) and x^ is known as 
the diffusion length, is a measure of the breadth of the CR precursor to the shock. It 
is important that computational resolution be fine enough to resolve the flow features on 
this scale. For our discussion we define, therefore, the resolution ratio of each simulation 
to be 

T = & < 2 - 12 > 

where Ax is the size of a computational zone. Frank, Jones & Ryu (1993) found that 
provided n r > 10 — 20, our code produced accurate and converged solutions. When 
ku 7^ kj_, td is not well defined, since k changes through a shock. To be specific, however, 
we will adopt the convention of using the value of k from equation [2.10] upstream of the 
shock to define in our discussions. 

Our numerical method solves equation [2.1] through the second order accurate, Eu- 
lerian finite difference, "Total Variation Diminishing" (TVD) scheme. The TVD scheme 
was originally designed for gasdynamics by Harten (1983) and extended to MHD by Brio 
and Wu (1988). In this paper, we use the version of the MHD-TVD code which was sub- 
sequently improved and fully tested by Ryu & Jones (1993). The pure MHD (S = 0) 
form of the equation is solved with the aid of an approximate MHD Riemann solver used 
to estimate the flux, F. CR source corrections, S, are added separately in a manner 
preserving second order accuracy. Shocks and other discontinuities are generally sharply 
resolved within a few zones. The CR energy equation is solved separately using a second 
order combined Lagrangian Crank-Nicholson and monotone remap scheme. This mixed 
scheme is efficient and provided us much better stability and accuracy than a straight 
Eulerian Crank-Nicholson approach. In the appendix we present the numerical method in 
more detail. A pure gasdynamical version of the code was tested extensively against both 
analytical steady state solutions and numerical time dependent models calculated with our 
well-tested PPM two-fluid code (Jones & Kang 1990), all with excellent agreement. The 
pure MHD code was tested using a nonlinear MHD Riemann solver against a variety of 
planar shock tube problems involving all three wave mode families (for details, see Ryu & 
Jones 1993). 

3. Results 



We will now summarize results from these simulations, beginning with some properties 
of flows that were continued until the state immediately downstream of the shock was no 
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longer evolving; i.e., a time asymptotic state. Then we will explore behaviors of CR-MHD 
shocks at intermediate times, while the post shock flow is still developing. There are also 
some notable features that continue to evolve downstream of the equilibrium shock. 



To describe the flows there are four parameters needed in addition to the CR "fluid" 
properties, 7 C and k. Thus, we parameterize our upstream state by the following dimen- 
sionless quantities, 



a 2 _ lg P g 



b 2 B 2 

Pr 

N = — , (3.1) 

Pa + Pc 



M f = 



9 t 1 c 

U X \ 



^ = tan -i(|i) 



In equation [3.1] all quantities refer to the upstream (incoming flow) state and u x 
is measured in the computed (piston) frame. In the following discussion wherever it is 
necessary to distinguish upstream and downstream conditions they will receive subscripts 
1 and 2 respectively. Of the remaining variables in equation [3.1] a is the gas sound 
(acoustic) speed, cj is the fast mode wave speed; namely, 



- = ,M. 0.2) 



A = hb 2 + a 2 + [{b 2 + a 2 ) 2 - 4b 2 a 2 cos 2 8 }^}, (3.3) 
J 2 

and 9 is termed the "obliquity" of the shock. Note that the fast mode Mach number, My, 
in our simulations refers to the piston rather than the shock. The equivalent shock Mach 
number is 

M fs = HM f . (3.4) 

I a x\ 

The relationship between Mj and Mj s depends on the compression through the shock, 
which is generally a function of the solution itself and a function of time as P c 2 changes. 
For most of the shocks described below, Mj < Mjg^/iMj. Our definition of (3 differs 
from that used frequently in the plasma physics literature and by WDV by a factor lg/2 
= 5 /6. For all the simulations described, the upstream flow is characterized by p = 1, 
u x = -1 and iV = 1/2. 

a) Steady CR-MHD Shocks: Comparison with Previous Calculations 



We start with a comparison between properties of a steady state model computed 
by WDV and the approximately equivalent time asymptotic flow found by us. Previously 
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(Frank, Jones, & Ryu 1993) we showed good agreement between our results and the special 
"cold gas" analytic solutions examined by Webb. These comparisons are useful not only 
as tests of numerical methods, but also as demonstrations that the analytic steady states 
can be reached from some other initial state. In that context, however, we do not address 
in this paper issues associated with multiple steady state solutions found by WDV for 
CR-MHD shocks (see Donohue et al. 1993 for such a discussion regarding CR-gasdynamic 
shocks) . 



Figure 1 illustrates the evolution to an apparent dynamical steady post shock state of 
a simulated flow with upstream parameters (3 = 1, 7 C = 4 /3, k = 10, and o = 60°. The 
resolution ratio was n r = 33, so that the solution is well converged. The simulation was 
designed to match as closely as possible a solution presented by WDV in which Mj s = 6.1. 
(Note that their definition of Cf included P c as a contributor to the sound speed. We have 
corrected for that in establishing an appropriate comparison.) The downstream partial 
pressures, P C 2, P g 2 and Pb2, from WDV are indicated by dashed lines. The final u s 
(and thus Mf s ) is hard to predict exactly, because it depends on P c 2, P g 2 and Pjs<i in 
nonlinear ways. In addition the WDV solution was read from a published figure (Figure 
12 in DWV), so we cannot match the two problems exactly. Nonetheless we find that all 
the partial pressures agree to within about 1% when normalized by the total momentum 
flux incident on the shock, P, defined as 

P = piv? 8 + P m + Pgl + Pd, (3-5) 

where to estimate u s in the final steady shock, we assume the mass flux jump conditions 
are exactly satisfied across the shock. Then the required u s in the frame of the upstream 
fluid can be written in terms of jump conditions across the shock in the computed frame 

as 

u s = -n — u x i = -Ki, ( 3 - 6 ) 

[Pi r - 1 

where 

r = (3.7) 

Pi 

is the compression ratio through the full shock structure. 



Figure 1 shows the flow at times t/t^ = 325, 856 and 1377. The shock transition has 
approached a steady state structure by t/t^ = 856. The rather large number of diffusion 
times needed for this shock to come into a dynamic equilibrium is consistent with previous 
gasdynamical CR shock calculations (e.g., Jones & Kang 1990; Kang & Jones 1991). 
Jones & Kang (1990) pointed out that the rate for CR pressure to increase in response to 
diffusive acceleration depends not only on tj, but also inversely on the CR adiabatic index. 
Physically that derives from the fact that pressure generated by non-relativistic particles 
scales as the particle momentum squared, while it increases only linearly with momentum 
for relativistic particles. Since 7 C = 4 /3, the pressure in this model is assumed to come 
entirely from relativistic particles. The time to reach equilibrium also depends on the total 
change in P c , of course, and so on N and Mj s . 
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As P c builds within the shock precursor, a postshock density spike forms and then is 
left behind after the shock approaches equilibrium. The origin of this feature in strongly 
modified gasdynamic CR shocks was explained by Drury (1987) and independently by 
Jones & Kang (1990). Briefly, it reflects an adiabatic compression in the precursor at 
times before an initially strong subshock weakens in response to substantial adiabatic 
heating in the precursor. The tangential magnetic field participates in the development 
of that feature in oblique MHD shocks, as can be seen in Figure 1. These features are 
examples of time dependent effects that occur in the evolution of CR modified MHD shock 
structures and of which the steady analytical models can give no hint. We will discuss 
these time-dependent features in more detail in section 3e. 



b) Steady CR-MHD Shocks: CR Acceleration Efficiency Maps 



We now explore the influence of tangential magnetic fields on the CR acceleration 
efficiency. All previous related discussions have pointed out that strong, oblique fields 
tend to inhibit acceleration efficiency, at least in shocks of modest strength. However, 
the physical reasons for this have not been well determined. Based on their study of 
magnetosonic shocks Jun et al. (1993) made the reasonable conjecture that efficiency was 
reduced because a tangential field inhibits compression by stiffening the equation of state 
for the fluid . But the physics of oblique CR-MHD shocks is complex and nonlinear, so this 
question needs further examination to be understood. For this discussion we will define 
the CR acceleration efficiency to be the ratio of downstream CR pressure to upstream total 
momentum flux, 

ec = 5*. (3.8) 

In their study, WDV used analytical solutions to the steady CR-MHD shock equations for 
Mj = 5 and 1 and O = 0°, 30° and 60° to conclude that the obliquity of the field was the 
important acceleration inhibiting factor. To extend that effort we have completed a large 
set of simulations designed to map out time asymptotic values of e c as a function of the 
upstream parameters. A compact, visual way to examine how e c depends on the magnetic 
field properties in the flow is to construct maps of e c on a plane defined by the in-flowing 
normal and tangential Alfvenic Mach numbers, 

M x = — = — = — ^— oc — , 3.8a 

b x b cos O cos 6*0 B x 

and 

M z = — = . = oc — , 3.86 

b z o sin O sin O B z 

where M g is the sonic (gas) Mach number. For each map j3 is held constant. Shocks with 
a range of shock fast mode Mach numbers, Mj s , are obtained by varying M x and M z . 
For this purpose it is convenient to write Mj in terms of (3 and the total Alfvenic Mach 
number, M , 

o 2M, 2 

Mf b - -, (3.9) 

J {1 + /3+ [(l + /3) 2 -4/3cos 2 o ]2} 
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where 

1 1 1 b 2 1 

It is also useful to note that 

M b = M sin 20 o , (3.11a) 
o = tan~ 1 (M x /M z ), (3.116) 

and 

M Z 2 = ^L ^ - 1 - 2/3, (3.11c) 

z Pg AT = 1/2 Pg p ' v ; 

where 

M = (M x + M|)3. (3.11d) 

Figure 2 presents such time-asymptotic efficiency maps for (3 = 20, 3, 0.8, and 0.2. 
Each map was made from 25 simulations which uniformly covered the (M x , M z ) plane to 
include 2 < My s < 16. All of the simulations have 7 C = 1.41 and n r > 20. The contours 
of CR acceleration efficiency, e c , are shown as solid lines and are labeled. The contours 
of fast mode shock Mach number, My s , are shown as dotted lines. They begin with the 
lowest contour set at My s = 2 in the lower left corner, increasing in steps of 2. Although 
the shocks in (3 < 1, My s ~ 2 flows would likely be influenced by finite "Alfven transport" 
effects {e.g., Jones 1993a), their inclusion here is useful in establishing trends within the 
CR-MHD solutions. 



Some preliminary explanation should make the physics illustrated in Figure 2 easier to 
grasp. Note first that moving along an arc of constant radius (M) from the vertical axis 
corresponds to sweeping through 9 from 0° to 90°. The total Alfvenic Mach number, M5, 
and the fast mode Mach number, My, are maximum on this arc at 9 = 45°. We note, 
as well, that for large M z contours of My asymptote to lines M x = My, while for large 
M x these contours asymptote to lines M z = ^/l + (3Mj-. Thus parallel and perpendicular 
shocks lie only at infinity in this plane. 

It is easiest to begin by examining the f3 = 20 map, since MHD effects are minimal. 
There one can see that both the e c and Mf s contours are almost symmetric about o = 45° 
and that they are almost "parallel". There is a weak tendency at large obliquity for 
e c to decrease along a Mach number contour as 6 Q increases. This is consistent with 
the conclusion reached by WDV. The effect is seen to be much stronger as one looks at 
successively smaller (3; i.e., as the field becomes dynamically more important. Consider, 
for example, the obliquity at which e c drops below 0.6 on the Mj s = 8 contour. As (3 
decreases from 20 to 0.2 the value of the upstream obliquity at which these contours cross 
drops from 9 ~ 63° to o ~ 47°. Thus, as the magnetic field becomes more dynamically 
significant shocks of constant Mj s and constant 9 become less efficient CR accelerators 
found by others. Relation [3.11c] provides a clue to what is happening. Lines of constant 
M z correspond to a constant ratio of in-flowing momentum flux to tangential magnetic 
pressure. When [3 is small Pb/P ^ M~ 2 . For [3 < 1, e c contours in Figure 2 for high 
obliquity (quasi perpendicular) shocks are closely spaced and roughly parallel to the M x 
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axis, and so parallel to constant Pg/P. Thus, it is the dynamical capacity of magnetic 
pressure that seems to be most significant. This is reasonable and consistent with previous 
conclusions, of course. The effect is important only for small to moderate Mj s , since as 
Mj s becomes large M z is restricted to values > y/1 + f3Mj 3> 1. Then from equation 
[3.11c] Pg/P must be small for any f3. This conclusion still does not discriminate among 
at least two possible physical reasons why the magnetic pressure plays this role. It could 
be, as suggested by Jun et al. (1993) just that the magnetic field reduces compression 
through the shock. Alternatively, it could reflect the need for the shock to do work on the 
field in compressing it, thus reducing the energy available to CR. 

We will describe an experiment in the following subsection that may help to resolve 
that issue. But, first we should comment on another aspect of the efficiency contours in 
Figure 2. In this case we focus on quasi parallel shocks, so that M x is relatively small. 
The efficiency of these shocks in the strong field (f3 < 1) regime seems to be higher at a 
given fast mode Mach number than when the field is weak. For instance, when 9 = 20° 
and (3 = 0.2, a Mj s = 4 shock leads to e c ~ 0.55. At the same angle this efficiency 
requires Mj s > 6 when (3 = 20. On the face of it this seems to imply that the strong field 
has actually enhanced e c . But that would be a misinterpretation of the results. When 
/3 <^ 1 Mj Ki i.e., the fast mode wave is much like an Alfven wave for quasi parallel 
propagation and not very compressive. The large value of B that makes cj ~ b so large 
is mostly in the normal component. Using M z as an indicator again, we can see that the 
dynamical impact of the tangential field cannot be very large in this regime. Thus, the 
fast mode Mach number underestimates the strength of the shock as it relates to diffusive 
acceleration. Indeed, examination of the structures of the CR-MHD shocks here reveals 
that the time asymptotic compression is large (r > 4) and comparable to that found in 
large (3 shocks having similar acoustic Mach number M g = M^j ' \f$. 

This last feature should emphasize that it would be wrong to expect that the existence 
of strong fields necessarily leads to a reduction in the efficiency of CR acceleration in 
astrophysical environments. That statement seems to apply only when the tangential 
magnetic pressure is strong. Thus, the orientation of the field is particularly important as 

(3<1. 

c) Steady Solutions: The Role of Compression and Magnetic Pressure 

Diffusive shock acceleration results from simultaneous compression and diffusion of 
CRs. Thus, the compression through the shock is itself obviously important to the efficiency 
of the process. But nonlinear effects complicate this simple statement. Even in CR- 
gasdynamical shocks the acceleration efficiency represents a balance in the competition 
between energy input into the CR and the gas. The incoming momentum flux does 
work on both the gas and the CRs, with the partitioning dependent on the required 
entropy production through the structure. In general a simple relationship between total 
compression and efficiency does not exist. When magnetic pressure is included the situation 
is even more complex, since compression of the magnetic field also requires work. 
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To explore the relative roles of the compression itself and the work done on the magnetic 
field in CR-MHD shocks, we have performed a simple set of numerical experiments. We 
carried simulations out to dynamical steady state postshock conditions for two values of 
the upstream obliquity, 9 Q = 30° and 60° and two values of (3 = 5, 1. For each of these 
combinations we computed solutions for a range of CR adiabatic index, 4 /3 < 7 C < 5 / 3 - 
By varying 7 C we change the compression through the shock, since 7 C < softens the 
equation of state in a manner analogous to the way that the magnetic pressure stiffens it 
compared to the gas alone. All of the models used My = 6. 

The results are illustrated in Figure 3. Figure 3a (upper left) shows the compression 
ratio, r, for the 9 = 60° models as a function of 7 C . As we expect from previous studies 
(e.g., Drury & Volk 1981; Achterberg et al. 1984; Kang & Jones 1990) smaller 7 C leads to 
substantially greater compression through the shock. Note, however, that the compression 
is independent of (3 within our levels of uncertainty. As seen in Figure 3b (upper right) 
the efficiency, e c , does depend on 7 C , and thus on the compression. Note, however, that e c 
actually decreases with increasing compression; i.e., as the mixed fluid develops a softer 
equation of state. That trend runs counter to the compression argument given by Jun et 
al. (1993) regarding the role of the magnetic field. But, e c also depends on (3. More to the 
point, for a fixed compression (7c), e c is reduced as the tangential magnetic field becomes 
stronger. The change in e c can be accounted for roughly in the work done on the field 
itself. That can be seen in Figure 3d, which shows the normalized downstream magnetic 
pressure, P^/P, for the same study. Two points are relevant. First, PbiIP varies in a 
sense opposite to e c as 7 C varies. Second, the difference in Pb2/P between the weak and 
strong field cases is comparable to that in e c for the same two cases. Thus, crudely, at 
least, there has been an exchange between these two quantities. But, the nonlinear nature 
of these processes would lead us to expect a more complex situation quantitatively. Indeed, 
for the smaller obliquity case, 9 = 30°, seen in Figure 3c, the efficiencies for different (3 
are the same to within our ~ 1% uncertainties. The magnetic pressure, PbiIP is smaller 
and has the same trend as in the other example. But, the low and high (3 values of PB2/P 
differ by as much as 3%. Clearly, the other energy channels (gas pressure and kinetic 
energy) are also involved. 

d) Time- dependent Effects: N anisotropic Diffusion 

From this point we will concentrate on evolutionary aspects of CR-MHD shocks as 
we can understand them through two-fluid models. For a variety of reasons it is likely 
that most astrophysically important CR shocks will not be of a steady character (see, e.g., 
Jones 1992). 

It is well known that the time asymptotic efficiency of two-fluid CR-gasdynamical 
shocks should be insensitive to the value or form of k (Drury & Volk 1981; Jones 1993b). 
On the other hand the detailed properties of k do influence the rate of shock evolution 
and the structure of the shock. We now briefly examine those issues for oblique CR-MHD 
shocks. Figure 4 presents properties of two simulations differing only in the properties 
of k. The other properties of the flows were: (3 = 1, Mj = 6, 7 C = 1.41, 9 = 60°. 
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Both shocks are shown at the same dynamical time, t, after the shock itself has reached 
an equilibrium structure. The dotted line represents a shock computed with constant 
k = 0.325, while the solid line presents a flow with an anisotropic (obliquity dependent) k 
with K|| = 1.0 = 10 x as defined in equation [2.10]. Upstream of the shock k = 0.325 in 
this model also. Two points are obvious from the figure. First, the final post shock states, 
as represented by p and P c , are the same in both simulations. This holds for the other 
state variables, as well. Second, the histories of the two shocks were different. One can 
trace the history by looking from left to right, since the shock moves away from the piston 
at the left boundary and to the first approximation the post shock fluid is at rest and not 
changing very fast (the next section discusses ways in which this is not quite true). The 
density spike identifies the approximate location of the shock when it reached dynamic 
equilibrium (see Figure 1). Thus, it is clear that this condition was achieved at an earlier 
time for the obliquity dependent diffusion model. That results from effective differences 
in diffusion times for the two models. For the constant diffusion model, = 0.325, which 
matches the nominal (upstream) value for the obliquity dependent diffusion model, using 
the convention established in §2. But, through the precursor of that shock k decreases 
from k = 0.325 to k = 0.1, so that effective values of tj and within the precursor and 
postshock flows are smaller as one proceeds downstream. That accounts in part for the 
sharper and higher post shock density spike seen in the variable diffusion model and the 
earlier achievement of a steady post shock state. Furthermore, the quicker evolution to a 
CR dominated post shock state with higher compression (due to the softer equation of state 
of the CR) explains why the position of the shock is not as far advanced for the obliquity 
dependent diffusion model. The shock speed decreases as the compression increases. 

The more rapid shock evolution in the obliquity dependent diffusion model is related 
to the point emphasized by Jokipii (1987). He argues from test particle theory that with 
such diffusion models acceleration of individual particles will be quicker in highly oblique 
shocks, raising the possibility that particles of higher energy might be achieved in supernova 
remnant shocks than expected with isotropic, Bohm diffusion. Since we have not solved 
the diffusion-advection equation here we cannot comment directly on that issue except to 
make one point. More rapid acceleration also leads more quickly to the development of 
nonlinear feedback from the CR to the fluid. In highly oblique shocks in low (3 plasmas, 
that would mean a prompter onset of work done by CRs on the tangential magnetic field. 
The net result of that, of course, is a reduction in the total work done by the flow on the 
CR population. Since for 7 C > 4 /3, P c is dominated by low energy particles, these two 
behaviors are not necessarily contradictory. Clearly, a close look at this issue is warranted, 
however. 

e) Time Dependent Effects: MHD Evolution 

In Figure 5 we present the time development of a Mj = 10, strong field ([3 = 1), CR 
modified oblique (9 = 60°) MHD shocked flow. The evolution of the shock was calculated 
with 7 C = 1.41 and nonisotropic diffusion using = 1 and Kj_ = 0.1. Three time frames 
are shown at t/t^ = 62, 123, and 185 respectively. The plots demonstrate that the coupled 
CR-MHD evolution of the shock and especially the post shock flow are complex. Basic 
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qualitative evolution of the fluid and the CRs are similar to CR-gasdynamical shocks, 
however, and they have been discussed extensively in previously cited literature. We 
concentrate here on those features that are MHD in origin. 

Through the shock the tangential magnetic field is enhanced. Mostly this is just 
compression, but in fast mode shocks there is an additional, Mach number dependent, en- 
hancement through induction generated by shear in the shock transition (e.g., Kantrowitz 
& Petschek 1966). That shear is evident in the u z jump seen through the shock. The net 
tangential magnetic field change through any steady plane compression feature is given by 
the simple expression 

B z2 ,Mj!-l, . . 

where M x refers to the condition upstream of the feature as measured in the rest frame of 
the feature. The factor on the right in brackets results from shear in the transition. The 
shock in Figure 5 corresponds to M x w 27, so rg r. But for weaker shocks, including 
slow shocks the shear induced term is more significant. Equation [3.12] also shows the 
well-known difference between fast mode and slow mode waves. Whereas the sign of the 
change in tangential field and density are the same for fast waves (Mj > r > 1), they are 
opposite for slow waves (Mj < 1). This difference results from a reversal in the direction of 
the shear, [u z ], relative to the normal velocity jump, [u x ] through the feature (Kantrowitz 
& Petschek 1966). Then the induction (see the magnetic terms in equation [2.1]) also 
changes sign. This information helps sort out the evolution of the postshock flows seen 
in strong CR-MHD oblique shocks such as that shown in Figure 5. In this regard it is 
also useful to identify from the MHD dispersion relation that through a f3 > 1 fast mode 
compression, (M x [u z ] ) / (M z [u x ] ) > 0. The opposite inequality applies for slow waves. 
For the flow described here M x /M z > 0, so the normal velocity gradient and the shear 
velocity gradient will have opposite signs in a fast compression, but the same sign in a 
slow compression. 

Although the shock itself in the simulations shown reaches a dynamical steady state 
that is not true for the post shock density (and magnetic field) spike. Such spikes result 
from initially strong subshocks as they adjust to new CR dominated equilibrium conditions. 
At the start this feature is a fast mode feature, since the tangential field and gas density 
increase together with similar form at early times. It develops, after all, from the fast 
mode shock. At early times one can see that the shear velocity component, u z that leads 
to inductive generation of field is negative behind the (fast mode) shock as expected from 
the relation just cited. But after the shock has reached an equilibrium, so that the density 
spike is left behind, a clear maximum in u z is present just to the left of the peak in the 
density spike. That signals the formation of a slow wave from the density peak. This 
interpretation comes from the fact that du x /dx < across almost the entire density spike, 
including the peak. Thus it is being slowly compressed. Hence the negative gradient in u z 
must represent a slow mode compression. This development is even more obvious by the 
third and last time shown in Figure 5. By this time, in fact, the (right facing) slow mode 
compression has substantially reduced B z on its downstream side, so that the density and 
magnetic spikes are clearly distinct. 
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Subsequently, the consequences of the slow wave are striking, in fact, as seen in Figure 
6. This figure shows at later times the evolution of the density spike. One sees that the 
magnetic spike is eventually confined to the forward (right) edge of the density spike. In 
fact the forward edge of the density spike has steepened into a slow mode shock. This is 
demonstrated by the sharp jump in u x at this point. The development of shocks in the 
transient density spike has not happened in any of the CR-gasdynamical shocks we have 
studied, so it seems to be strictly the result of the CR-MHD nature of this calculation. 
Although the density spike results from artificial nonequilibrium initial conditions, it may 
still be of much more than academic interest. Many astrophysical shocks will likely have 
transient properties. We have seen such density spikes in our simulations of SNR (Jones 
& Kang 1992) and in 2D simulations of supersonic dense clumps (Jones & Kang 1993; 
Jones, Kang & Tregillis 1993). The velocity jump at the slow mode shock in the spike 
is moderately large. It actually does lead to a small amount of CR particle acceleration 
downstream of the main shock. That can be seen through the small "bump" in the P c 
curve at the last time shown in Figure 6. In addition, the shear and associated tangential 
current in the postshock flow resembles in some ways those inside the main shock. There 
they are accompanied by drift acceleration. For calculations in a planar geometry this does 
not appear, because the net energy flux vanishes (Jokipii 1982; Webb, Axford & Terasawa 
1983). However, once that symmetry is broken, the drifts could play important roles just 
as they apparently due in perturbed oblique CR-MHD shocks (Decker 1988; Zank et al. 
1990). 

4. Summary and Conclusions 

Utilizing a new numerical method built on the two-fluid model for diffusive shock 
acceleration we have conducted a study of CR modified oblique MHD plane shocks. The 
results of our numerical experiments are as follows: 

1. For a fixed fast mode Mach number the particle acceleration efficiency, e c , of time 
asymptotic CR modified quasi perpendicular MHD shocks decreases as the upstream 
obliquity, O , of the magnetic field increases. The dependence of e c on 9 is greater for 
tangential fields whose amplification through the shock transition absorbs a significant 
portion of the incident momentum flux of the CR-MHD fluid. This limits the impact of 
the field on acceleration efficiency to moderately strong, quasi perpendicular shocks in 
strong field (low f3) plasmas. These results confirm early analytical models and recent 
numerical models of perpendicular, magnetosonic shocks. Our experiments, however, 
map the relationships between e c and the MHD and CR properties more fully than 
previous studies and demonstrate that acceleration inhibition probably results from 
work done on the magnetic field rather than changes in fluid compression brought 
about by the presence of the field, for example. 

2 For quasi parallel shocks, our results show that e c is higher at a given fast mode Mach 
number when (3 < 1 than when (3 > 1. While this would appear to suggest that strong 
fields enhance acceleration efficiency in quasi parallel shocks it actually results because 
in strong fields the quasi parallel fast wave speed is the Alfven speed. The larger gas 
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sonic Mach number is a better measure of the strength and therefore the compression 
of (3 < 1 quasi parallel shocks. Thus strong fields do not necessarily imply a reduction 
in e c . It is the tangential magnetic pressure that is important. 

3. The possibility of an anisotropic form for the CR diffusion coefficient, k, does not affect 
the time asymptotic post shock state achieved by the CR modified MHD shocks. How- 
ever, the history of CR-MHD shocks can be considerably modified by different forms 
of k. In particular we have examined the effect of an obliquity dependent diffusion 
coefficient k{6) on CR-MHD shock evolution. In these models the compression of the 
tangential field through the shock reduces k in the downstream regions and reduces 
the time necessary to achieve a steady state. That is consistent with earlier sugges- 
tions that anisotropic diffusion in quasi perpendicular shocks might allow individual 
cosmic ray particles to reach higher energies in astrophysical shocks having finite life 
times. The fact that strong field, quasi perpendicular shocks may at the same time 
become less efficient accelerators as the result of nonlinear effects needs to be studied 
using more complete models of diffusive acceleration. The two outcomes need not be 
contradictory, however, since the efficiency measures total energy exchange rather than 
the distribution of energies. 

4. Time dependent CR-MHD effects may become significant downstream of the shock, 
even once the shock itself has approached a dynamical equilibrium. In particular if 
CR modification of the shock is sufficient to generate a postshock density spike, MHD 
influences in that spike can cause substantial evolution in the spike. We found that 
a strong post shock shear develops in response to magnetic stresses in strong, low (3 
oblique shocks. In more realistic flows where planar symmetry constraints are broken 
this shear layer may yield additional CR acceleration through the production of drift 
currents. Further, a slow mode shock can form on the near edge of the post shock 
spike, possibly enhancing particle acceleration from the full structure. 
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Appendix: Finite Difference Algorithms 



To solve the MHD equations we have employed the Total Variation Diminishing (TVD) 
method introduced by Harten (1983). The TVD method solves the vector equation [2.1] 
without the source vector S in conservative, finite difference form: 

( 9 r i r = 9 ?+^(/» +i «)-?»,(«). (a.i) 

As usual, i is the spatial zone index; integers correspond to zone averages (centered values). 
The n index refers to time step. fV^i(q) is a 4 point numerical flux function that depends 

on the vector q and the original flux F(q). The TVD prescription of the numerical flux 
function f n ,i(q) utilizes an approximate MHD Riemann solver (Roe 1981; Harten 1983; 

Brio & Wu 1988; Ryu & Jones 1993) to find estimates of the state variables at the zone 
boundaries. The method is an explicit, second order, Eulerian finite difference scheme. The 
TVD scheme also applies minimal numerical viscosity to the computed solution resulting 
in strong shocks that are represented with only 2 to 3 zones. Contact and rotational 
discontinuities are generally contained within <10 zones. For detailed descriptions and 
tests for the TVD method applied to the MHD problem see Ryu & Jones (1993). 

CR dynamical feedback on the conducting fluid is handled in a two step process 
designed to preserve second order accuracy. Gradients in P c accelerate and do work on the 
gas during a time step, At. To reflect that in the fluxes, f^_i(q), and thus improve the 

performance of the MHD Riemann solver, input Riemann solver momentum and energy 
densities are corrected to time centered values by a simple estimate of the force and work 
produced by the gradient in P c , 

Q = Q? + ^Sf. (A.2) 

Then, after the TVD step represented in equation [A.I] with q, the state vector q is 
updated for the effects of the source terms, 5; namely, 

q? +1 = (q? +1 T+AtSf, (A.3) 

where q* represents the solution to equation [A.I]. 

The CR energy is then updated through a two step procedure, which is also second 
order accurate. In the first step we solve the Lagrangian version of equation [2.5], 

dE c 8E C dE c d ( dE c \ ^ du n .. „. 

exactly as we do in the PPM two-fluid code, using a Crank-Nicholson scheme explained in 
Jones & Kang (1990). This is followed by a "remap" step, 

Kt l = (*#V - «?^G£ + i " 9U 2 (E*)), (A.5) 
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where 

g» ±h (E* c ) = (E^y ± ^{E n c /{1 ± (A.6) 

E* is the solution to equation [A. 5] and (E c )* is the monotone slope of E*, (Van Leer 
1974). In equations [A. 4 - A.6] the x subscript has been dropped from u x to simplify 
notation. This operator split method of solving equation [2.5] proved to offer the best 
stability and accuracy among several approaches we tried in numerical experiments during 
code development. 
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Figure Captions 

Fig. 1. — Comparison between the analytical post-shock state (dashed line) and the results 
of a time-dependent simulation with approximately the same upstream conditions. 
Both models have Mf s = 6.1, 9 = 60°, 7 C = 4 /3 and (3 = 1.0. The time-dependent 
simulation used constant k = 10. Shown are the density, gas pressure, magnetic 
pressure, normal velocity, tangential velocity and CR pressure. The simulation is 
shown at t/t^ = 325,856 , and 1377, where the diffusion time, tj is defined in the 
text. 

Fig. 2. — Maps of CR acceleration efficiency, e c , as a function of normal and tangential 
Alfvenic piston Mach number (M x , M z ) for four different values of (3. Solid lines are 
contours of e c , dashed lines are contours of fast mode shock Mach number, Mj s . e c 
contour levels are marked. In all the maps contours of Mf s are separated by steps of 
2, with Mf s = 2 always the first contour seen in the lower left corner of map, and 
Mj s = 16 in the upper right. 9 increase towards the x axis. 

Fig. 3. — Results from a series of Mf s = 6 runs with (3 = 5 (solid line) and (3 = 1 
(dashed line) for a range in the CR adiabatic index, 7 C . Upper Left: Compression 
ratio, r = p2/ ' p\, as function of 7 C for upstream 9 = 60°. (3 = 5 (solid) and (3 = 1 
(dashed). Upper Right: CR acceleration efficiency e c as function of 7 C for upstream 
9 = 60°. Lower Left: CR acceleration efficiency e c as function of 7 C for upstream 
9 = 30°. Lower Right: Normalized downstream magnetic pressure, as function 
of 7 C for upstream 9 Q = 60°. Note that the (3 = 1 case is always more sensitive to 7 C 
and, therefore, to compression. 

Fig. 4. — Comparison of two (3 = 1, Mt = 6, 7 C = 1.41 runs and different models of the 
CR diffusion coefficient. Dashed line shows a constant n = 0.325 model. The solid line 
shows an obliquity dependent model (see text) Both represent the shock structures at 
t/t d = 431. 

Fig. 5.— The time evolution of a model with M f = 10, 9 Q = 60°, 7 C = 1.41, P go = 0.032 
and P = 1.0. k is obliquity dependent as in Figure 4. Structure is shown at times 
t/td = 62, 123 and 185. Shown are the density, gas pressure, magnetic pressure, 
normal velocity, tangential velocity, tangential field, obliquity, diffusion coefficient 
and CR pressure. 

Fig. 6. — Same as Figure 5 at times t/t^ = 308,615 and 923. Note that here we have 
changed scale in x to focus on the evolution of the post shock density spike. 



